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Abstract 

Weakly Interacting Massive Particles (WIMPs) are one of the leading candidates for 
Dark Matter. We develop a model-independent method for determining the mass m x of the 
WIMP by using data (i.e., measured recoil energies) of direct detection experiments. Our 
method is independent of the as yet unknown WIMP density near the Earth, of the form of 
the WIMP velocity distribution, as well as of the WIMP-nucleus cross section. However, 
it requires positive signals from at least two detectors with different target nuclei. In a 
background-free environment, m x ~ 50 GeV could in principle be determined with an error 
of ~ 35% with only 2 x 50 events; in practice upper and lower limits on the recoil energy 
of signal events, imposed to reduce backgrounds, can increase the error. The method also 
loses precision if m x significantly exceeds the mass of the heaviest target nucleus used. 



1 Introduction 



First indications for the existence of Dark Matter were already found in the 1930s [1]. By now 
there is strong evidence [l]-[5] to believe that a large fraction (more than 80%) of all matter in the 
Universe is dark (i.e., interacts at most very weakly with electromagnetic radiation and ordinary 
matter). The dominant component of this cosmological Dark Matter must be due to some yet 
to be discovered, non-baryonic particles. Weakly Interacting Massive Particles (WIMPs) x are 
one of the leading candidates for Dark Matter. WIMPs are stable particles which arise in several 
extensions of the Standard Model of electroweak interactions. Typically they are presumed to 
have masses between 10 GeV and a few TeV and interact with ordinary matter only weakly (for 
reviews, see [6]). 

Currently, the most promising method to detect many different WIMP candidates is the 
direct detection of the recoil energy deposited in a low-background laboratory detector by elastic 
scattering of ambient WIMPs on the target nuclei [7, 8]. The recoil energy spectrum can be 
calculated from an integral over the one-dimensional WIMP velocity distribution fi(v), where v 
is the absolute value of the WIMP velocity in the laboratory frame. If this function is known, the 
WIMP mass m x can be obtained from a one-parameter fit to the normalized recoil spectrum [9], 
once a positive WIMP signal has been found. However, this introduces a systematic uncertainty 
which is difficult to control. We remind the reader that iV— body simulations of the spatial 
distribution of Cold Dark Matter (which includes WIMPs) seem to be at odds with observations 
at least in the central region of galaxies [10]; this may have ramifications for j\ as well. On 
the theory side, several modifications of the standard "shifted Maxwellian" distribution [6] have 
been suggested, ranging from co- or counter-rotating halos [11] to scenarios where the three- 
dimensional WIMP velocity distribution gets large contributions from discrete "streams" with 
(nearly) fixed velocity [12, 13]. 

The goal of our work is to develop model-independent methods which allow to determine 
m x directly from (future) experimental data. This builds on our earlier work [14], where we 
showed how to determine (moments of) fi(v) from the recoil spectrum in direct WIMP detection 
experiments. In this earlier analysis m x had been the only input required; one does not need to 
know the WIMP-nucleus scattering cross section, nor the local WIMP density. The fact that 
this method will give a result for fi for any assumed value of m x already tells us that one will 
need at least two different experiments, with different target nuclei, to model-independently 
determine the WIMP mass from direct detection experiments. To do so, one simply requires 
that the values of a given moment of j\ determined by both experiments agree. This leads 
to a simple expression for m x , which can easily be solved analytically; note that each moment 
can be used. An additional expression for m x can be derived under the assumption that the 
ratio of scattering cross sections on protons and neutrons is known. This is true e.g., for spin- 
independent scattering of a supersymmetric neutralino, which is the perhaps best motivated 
WIMP candidate [6]; in this case the spin-independent cross section for scattering on a proton 
is almost the same as that for scattering on a neutron. Not surprisingly, the best result obtains 
by combining measurements of several moments with that derived from the assumption about 
the ratio of cross sections. 

The remainder of this article is organized as follows. In Sec. 2 we review briefly the methods 
for estimating the moments of the velocity distribution function, paying special attention to 
experimentally imposed limits on the range of allowed recoil energies. In Sec. 3 we will present 
the formalism for determining the WIMP mass. Numerical results, based on Monte Carlo 
simulations of future experiments, will be presented in Sec. 4. We conclude in Sec. 5. Some 
technical details of our calculation will be given in an Appendix. 



2 Determining the moments of the velocity distribution 
of WIMPs 



In this section we review briefly the method for estimating the moments of the one-dimensional 
velocity distribution function fi of WIMPs from the elastic WIMP-nucleus scattering data. 
We first discuss the formalism, and then describe how it can be implemented directly using 
(simulated, or future real) data from direct WIMP search experiments. 



2.1 Formalism 



Our analysis starts from the basic expression for the differential rate for elastic WIMP-nucleus 
scattering [6]: 
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Here R is the direct detection event rate, i.e., the number of events per unit time and unit mass 
of detector material, Q is the energy deposited in the detector, F(Q) is the elastic nuclear form 
factor, and v is the absolute value of the WIMP velocity in the laboratory frame. The constant 
coefficient A is defined as 
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where po is the WIMP density near the Earth and ctq is the total cross section ignoring the form 
factor suppression. The reduced mass m T ^ is defined as 



m T N = 



(3) 



where m x is the WIMP mass and mN that of the target nucleus. Finally, v m i n is the minimal 
incoming velocity of incident WIMPs that can deposit the energy Q in the detector: 

v min = a^Q, (4) 
where we define 



a 



(5) 



Eq.(l) can be solved for fi [14]: 
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where the normalization constant M is given by 
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Note that, first, because fi(v) in Eq.(6) is the normalized velocity distribution, the normalization 
constant M here is independent of the constant coefficient A defined in Eq.(2). Second, the 
integral here goes over the entire physically allowed range of recoil energies, starting at Q — 0. 
The upper limit of the integral has been written as oo. However, it is usually assumed that the 



WIMP flux on Earth is negligible at velocities exceeding the escape velocity t> CS c — 700 km/s. 
This leads to a kinematic maximum of the recoil energy, 
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where a has been given in Eq.(5). Eq.(7) then implies 
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Using Eq.(6), the moments of f\ can be expressed as [14] 

rv esc ( a n+1 \ 
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which holds for all n > 0. Here the integral I n is given by 
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In this notation, Eq.(7) can be re-written as J\f = 2/(aI ). 

The results in Eqs.(6) and (10) depend on the WIMP mass m x only through the coefficient a 
defined in Eq.(5). Evidently any (assumed) value of m x will lead to a well-defined, normalized 
distribution function f\ when used in Eq.(6). Hence m x can be extracted from a single recoil 
spectrum only if one makes some assumptions about the velocity distribution f\{v). 

A model-independent determination of m x thus requires that at least two different recoil 
spectra, with two different target nuclei, have been measured. As we will show in detail in the 
next section, m x can then be obtained from the requirement that these two spectra lead to the 
same moments of f\. 

Before coming to that, we have to incorporate the effects of a finite energy acceptance of 
the detector. Any real detector will have a certain threshold energy Qthre below which it cannot 
register events. Off-line one may need to impose a cut Q > Q m i n > Qthrc in order to suppress 
(instrumental or physical) backgrounds. Similarly, background rejection may require a maximum 
energy cut, Q < Q max ^ Qmax.km- in fact, we will see below that, at least for smallish data 
samples, such a cut might even be beneficial for the determination of m x . We therefore now give 
expressions for the case that only data with Qmin < Q < Qnmx are used in the analysis. 

To that end, we introduce generalized moments of f\: 
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where we have introduced the short-hand notation 
dR 

Q=Qi 



r(Qi) 



dQ 



i = l,2, 



and 



/n(Qi,Q 2 )= f Q2 Q (n - 1)/2 

JQi 




[F 2 (g) \dQJ\ 



dQ. 



(12) 



(13) 



(14) 



In order to arrive at the final expression in Eq.(12) we used Eq.(6) and integrated by parts. The 
Qi in Eqs.(12)-(14) are related to the original integration limits Vi appearing on the left-hand 
side of Eq.(12) via Eq.(4), i.e., 

g, = 4' i = l, 2. (15) 

Of course, Eq.(12) reduces to Eq.(lO) in the limit v\ — » 0, v 2 — > f C s C ; note, however, that Eq.(12) 
is also applicable for the case n — — 1, where the last term on the right-hand side vanishes. The 
same restriction on the WIMP velocity can also be introduced in the normalization constant M 
of Eq.(7), in which case fi(v) dv is normalized to unity. Formally this can be treated using 
Eqs.(12)-(15) by demanding (v°)(v 1 ,v 2 ) = 1. 



2.2 Experimental implementation 

In order to directly use our results for fi(v) and for its moments (v n ) given in Eqs.(6), (10) 
and (12), one needs a functional form for the recoil spectrum dR/dQ. In practice this results 
usually from a fit to experimental data. However, data fitting can re-introduce some model 
dependence and makes the error analysis more complicated. Hence, expressions that allow to 
reconstruct f\ (v) and its moments directly from the data have been developed [14] . We started 
by considering experimental data described by 

Qn ~ b T < Qn, <Qn+ b T, * = h % ' ' ' , N n , U = 1, 2, • • • , B. (16) 

Here the total energy range has been divided into B bins with central points Q n and widths 
b n . In each bin, N n events will be recorded. Note that we assume that the sample to be 
analyzed only contains signal events, i.e., is free of background, and ignore the uncertainty on 
the measurement of the recoil energy Q. Active background suppression techniques [15] should 
make the former possible. The energy resolution of most existing detectors is so good that its 
error will be negligible compared to the statistical uncertainty for the foreseeable future. 

Since the recoil spectrum dR/dQ is expected to be approximately exponential, we used the 
following ansatz for the spectrum in the n— th bin [14]: 

I = ( ^\ ~ r P MQ-Qn) = r p MQ-Q s ,„) (-i 7 ) 

Here r n is the standard estimator for dR/dQ at Q = Q n , 
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r n is the value of the recoil spectrum at the point Q = Q n , 
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and k n is the logarithmic slope of the recoil spectrum in the n— th bin. It can be computed 
numerically from the average Q— value in the n— th bin: 

knbn \ 1 



0_ .|.= f coth ^ (20) 



where 
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Finally, Q s ^ n is the shifted point at which the leading systematic error due to the ansatz in 
Eq.(17) is minimal [14], 
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Note that Q s ,n differs from the central point of the n— th bin, Q n . 

Using Eqs.(12) and (13) with Q\ = Qmin and Q2 = Qmax, the generalized n— th moment of 
the velocity distribution function can be written as 
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where t> (Q) = a\[Q. Here we have implicitly assumed that Q max is so large that terms oc r(Q max ) 
are negligible. We will see later that this is not necessarily true for n > 1, since these moments 
receive sizable contributions from large recoil energies [14]. Nevertheless we will show in Sec. 3 
that even in that case, Eq.(23) can still be used for determining m x . From the ansatz Eq.(17), 
the counting rate at Q m m can be expressed as 



r(Qmin) = ri e kl ^- Q ^ . 
The integral I n (Qmm, Qm&x) defined in Eq.(14) can be estimated through the sum: 
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where the sum runs over all events in the data set that satisfy Q a G [Qmm? Qm&x]- 
Since all /„ are determined from the same data, they are correlated, with [14] 
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where the sum again runs over all events with recoil energy between Qmm and Qmax- 
On the other hand, the statistical error of r(Q m - m ) can be obtained from Eq.(24) as 
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The error on r\ follows directly from its definition in Eq.(18): 
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The error on the logarithmic slope k\ can be computed from Eq.(20): 
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Finally, the correlation between the errors on r(Q mm ), which is calculated entirely from the 
events in the first bin, and on /„ is given by [14] 



COv(r(<2min),4) 

'"(Qmin) -^n(Qmin, Qmin + 6l) 



X 



V(n) 



+ 



x 



In+2 (Q mm? Ynun 
-^n(Qmiri) Qmin + 6l) 



coth 



'61 fci 



a 2 (fc0 



(31) 



note that the integrals U in Eq.(31) only extend over the first bin, which ends at Q — Q ni \ n + b\. 



3 Determining the WIMP mass 

We are now ready to describe methods to extract the WIMP mass m x from direct detection 
data. Recall that m x is an input in the determination of f\ and its moments as described 
in the previous section. In particular, m x appears in the factor a n on the right-hand side of 
Eq.(23) describing the experimental estimate of the moments of f±. A truly model-independent 
determination of m x from these data will therefore only be possible by requiring that two (or 
more) experiments, using different target nuclei, lead to the same result for j\. 

Here we focus on the moments of the distribution function, rather than the function itself, 
since non-trivial information about the former can already be obtained with O(20) events [14]. 
We will also show how m x can be estimated from the knowledge of the integral Io appearing in 
the normalization of f\, if the ratio of WIMP scattering cross sections on protons and neutrons 
is known. For greater clarity, we will first discuss the idealized situation where all signal events, 
irrespective of their recoil energy Q a , are included. In the second subsection we will introduce 
upper and lower limits on Q a . A third subsection is devoted to a discussion how the different 
estimators for the WIMP mass can be combined using a x 2 fit. 



3.1 Without cuts on the recoil energy 

If no cuts on the recoil energy need to be applied, we can use the original moments, or equiva- 
lently, the generalized moments with v\ — > and v 2 — > v csc ; recall that the latter is the same as 
allowing t> 2 — * 00, since by assumption fi(v) = for v > v esc . Eq.(23) then simplifies to 

(v») = a n (n + 1) (|) , (32) 

where I n and I can be estimated as in Eq.(25) with Q m i n — > and Q max — > 00 (or, equivalently, 

Qmax ^ Qmax,kin)' 

Suppose X and Y are two target nuclei. We denote their masses by mx, my, and their form 
factors as Fx(Q), Fy(Q)- Similarly, we define axy as in Eq.(5), with mjy — > mx,Y- Eq.(32) 
then implies 

(33) 



Note that the form factors in the estimates of I ny x and I n y using Eq.(25) are different. Eq.(33) 
can be solved for m x using Eqs.(5) and (3): 

m x = , (34) 

ll n - ^mx/my 

where we have defined 

n n ^^=( I f^- I ^) 1/n , n^O, -1. (35) 

Using standard Gaussian error propagation, the statistical error on m x estimated from 
Eq.(34) is 
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where a 2 (I n ,x) = cov(I n ,x, In,x) and cov(/o,x, In,x) and so on can be estimated from Eq.(26). 

A second method for determining m x starts directly from Eq.(l), plus an assumption about 
the relative strength for WIMP scattering on protons p and neutrons n. The simplest such 
assumption is that the scattering cross section is the same for both nucleons. This is in fact 
an excellent approximation for the spin-independent contribution to the cross section of super- 
symmetric neutralinos [6], and for all WIMPs which interact primarily through Higgs exchange. 
Writing the "pointlike" cross section a of Eq. (2) as 

^o= (^m r 2 )N A 2 |/ P | 2 , (37) 

where f p is the effective XXPP 4-point coupling, m rj N is the reduced mass defined in Eq.(3) and 
A is the number of nucleons in the nucleus, we have from Eqs.(l), (2) and the first Eq.(12): 

r(Qmin) = ^(^Ayj 2 F 2 (Q mi A( v ^)( v (Q min ) :V ^ c ) . (38) 

Using the second Eq.(12) we see that the counting rate at Qmm m fact drops out, and we are 
left with 

i = M fjahiL.) , (39) 

where we have assumed Q m i n = 0. Recall that the rate dR/dQ is defined as rate per unit mass 
and observation time interval, i.e., we need to divide the actual event rate by the exposure 
£ = Mr, where M is the (fiducial) mass of the detector and r the observation time. In our 
previous discussion this factor always dropped out in the end, due to the appearance of the 
normalization N '. This is true also for the right-hand side of Eq.(38), but not for the 1/S factor 
on the left-hand side of this equation; hence a factor £ appears in the numerator of Eq.(39). 



Note that £ is dimensionless in natural units. On the other hand, the unknown factor po|/p| 2 
appearing in Eq.(39) will cancel out when we use this identity for two different targets X and 
Y, leading to the final result 
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Recall that, even though the derivation started from the expression for the counting rate at 
Q m in, this rate actually dropped out when going from Eq.(38) to Eq.(39). The final expression 
only depends on the quantity R a , which is estimated from all the events in both samples using 
Eq.(25). The error on m x computed from Eq.(40) is therefore comparable to that for m x derived 
from a moment of f\. Still keeping Q m i n = 0, we have 
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Ultimately the 7Z n , 7Z a and the errors in Eqs.(36) and (42) should be estimated from the data 
directly. In the meantime it is instructive to note that the final expressions for the statistical 
errors of our estimators for m x decompose into two factors. The expectation value of the first 
factor does not depend on it can be computed entirely from the masses of the involved 
particles: 
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Here we have made use of the identities 7Z n = a Y /a x in Eq.(35) and 
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which hold for the expectation values of these quantities as can be seen from Eqs.(34) and (40). 
Remarkably, the expectation values of the factors in front of the expressions in square brackets 
are in fact the same in Eqs.(36) and (42), apart from the factor l/\n\ appearing in the former 
equation. On the other hand, the expressions inside these square parentheses do depend on 
fi{v), as well as on the involved masses and form factors. 

It is nevertheless instructive to study the behavior of the factor k for different target nuclei X 
and Y . This largely determines what choice of targets is optimal, i.e., minimizes the statistical 
errors of our estimators of m x . It is easy to see that the ratio K,/m x , which will appear in the 
relative error a(m x )/m x , only depends on the dimensionless ratios m x /m x and m Y /m x . In 
Fig. 1 we therefore show contours of K,/m x in the plane spanned by these two ratios, using the 
ordering m Y > m x . We note first of all that k diverges as m Y — > m x . This is clear from the 
fact that the denominator in the final expression in Eq.(43) vanishes in this limit, while the 
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Figure 1: The five thick lines show contours of constant n/m x in the plane spanned by mx/fn x 
and my/m x , where we have taken my > mx without loss of generality. Here k, has been defined 
in Eq.(43); recall that the final relative error on m x is directly proportional to k/tti x . The lower 
thin line indicates the end of the physical region, my = mx, whereas the upper thin line shows 
my = (76/28)mx, corresponding to Silicon and Germanium targets. 

numerator is finite. Physically this simply means that performing two scattering experiments 
with the same target will not allow one to determine m x . 

Slightly less trivially, we also see that n/m x , and hence the relative error on m x , becomes 
very large if both target nuclei are either much heavier or both much lighter than the WIMP. 
This can be understood from the fact that m x only enters via the reduced mass defined in 
Eq.(3).* If m x 3> m^, m T ^ — > becomes completely independent of the WIMP mass. In 
the opposite extreme, m x <C m N , m r N — > m x becomes independent of the mass of the target 
nucleus, i.e., one is effectively back in the situation where both experiments are performed with 
the same target. 

These considerations favor chosing the two targets to be as different as possible. However, 
there are limits to this. On the one hand, taking a very light target nucleus will lead to a low 
event rate for this experiment, and hence very large statistical errors. Since the errors of both 
experiments are added in quadrature inside the expressions in square parenthesis in Eqs.(36) 
and (42), this would lead to a large overall error on m x . In fact, if the total event number is 

*Eq.(38) seems to have additional m x dependence. However, this comes from the factor n x = j o / m x' which 
drops out when the ratio of two targets is considered. 



held fixed, the final error on m x will be minimal if both experiments contain approximately the 
same number of events. 

On the other hand, taking a very heavy target nucleus also leads to problems. Heavy nuclei 
are large, which means they have quite soft form factors. For example, the Woods-Saxon 
prediction [16] for the form factor for 136 Xe, which is the target in some existing experiments 
[15], has a zero at Q ~ 95 keV. For our default choice v esc = 700 km/s, this is below Q max ,kin of 
Eq. (8) for all m x > 45 GeV. This is a serious problem, since the form factor, and hence the event 
rate, remain (very) small even beyond this zero. Experiments with 136 Xe therefore effectively 
impose a cut Q max < 95 GeV.* 

From these considerations it seems that chosing X = Si, Y = Ge might be a good compromise. 
This corresponds to points on the upper, thin solid straight line in Fig. 1. We see that in this 
case k/itl x > 4 for all values of m x . Not surprisingly, k/itl x reaches its minimum when m x lies in 
between the masses of the two target nuclei, i.e for the case at hand, for m x ~ 50 GeV. Note also 
that K,/m x is unfortunately always well above unity. One will thus have to get fairly accurate 
estimates of the relevant integrals 7j if one wishes to determine m x to better than a factor of 
two. 

Before discussing the statistical uncertainty in more detail, we proceed to include non-trivial 
upper and lower cuts on the recoil energy in our analysis. 

3.2 Incorporating cuts on the recoil energy 

As noted earlier, real experiments will probably have to impose both lower and upper limits 
on the recoil energy, partly for instrumental reasons, and partly to suppress backgrounds. This 
led us to introduce generalized moments of f\ in Eq.(12). These moments determined by two 
different detectors will still be the same, if either the integration limits are identical, v^x = Vi,Y, 
% — 1, 2, or the cuts are so loose that the contribution from WIMPs with v < v± or v > v 2 is 
negligible; a combination of these two possibilities can also occur, e.g., small but not necessarily 
equal values of Vi t x and v i ; y, and v 2 ,x — v 2 ,y significantly below v esc . 

The crucial observation that forms the basis for our analysis of m x as estimated from the 
generalized moments is that in fact all three terms in the final expression in Eq.(12) will agree 
separately between two targets, as long as both integration limits coincide. This can be shown 
by replacing the r(Qi) in Eq.(12) via Eq.(l) and using Eq.(15). In principle we could therefore 
simply equate the last (integral) terms between the two targets. This would lead to expressions 
very similar to those derived in the previous subsection, the only difference being that all integrals 
or sums over recoil energies would run over limited ranges. 

The problem with this approach is that the velocities Vi appearing in Eq.(12) are not directly 
observable. The recoil energies are. However, the values of Qmin and Qmax would need to be 
different for the two targets, if they are to correspond to the same values of v\ and v 2 . Worse, 
Eq.(15) shows that the ratios Q m m,x/Qmm,Y and Q m a,x,x/Qm&x,Y that one has to chose in order 
to achieve v it x = v%,y depend on m x . We are thus faced with the situation that, in the presence 
of significant cuts on the recoil energy, the quantity we wish to determine is needed as an input 
at an early stage of the analysis! 

This statement is, strictly speaking, true; we see no way around this, if significant cuts on 
the recoil energy are indeed needed. However, we can alleviate the situation. To begin with, we 
noticed that the systematic error one introduces by simply setting Q m in,x = Qmm,Y is reduced 
greatly if one keeps the sum of the first and third terms in Eq.(12). Moreover, we will show 

^Note also that while the integrals /„ remain finite where the form factor vanishes, the estimates for the errors 
will diverge there, due to the factor 1/F 4 appearing in Eq.(26). 



that, to some extent at least, one can do the matching of two Qmax values for the two targets 
directly from the data. In principle we could reduce the systematic error associated with the 
choice of Q max even more by also including the second term in Eq.(12). However, with limited 
statistics the error on r(Q max ) will always be very large, so that keeping this term will not help 
significantly 

Our determination of m x from generalized moments of fi is thus based on equating the sum 
of the first and third terms in Eq.(12) for two different targets. The resulting expression for m x 
can still be cast in the form of Eq.(34), but lZ n is now given by 



ZQ^TrxiQrnin^/FUQrn^x) + (n + l)I n ,X 



^Qmm,X r x(Qmm,x)/Fx(Q min! x) + h. 



X 



1/n 



(45) 



Here n 7^ and r^x,Y)(Qmin,(x,Y)) Te ^ eT to the counting rate for detectors X and Y at the 
respective lowest recoil energies included in the analysis. They can be estimated from Eq.(24); 
of course, the values of r%, k\, Q St ±, and bi will in general differ for the two targets. Note that 
the integrals I n , Jo are now to be evaluated with Q m i n as lower and Q max as upper limits, as in 
Eq.(25). 
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Figure 2: Ratio of reconstructed to input WIMP mass for imperfect Q m i n matching, for two 
different values of Qmin- We have taken Si and Ge targets and assumed infinite statistics and 
(effectively) no upper cut on the recoil energy. The solid (dashed) curves have been computed 
from Eqs.(34) and (45) with n = 1 including (neglecting) all terms oc r(Q min ). See the text for 
further details. 



Before discussing the statistical uncertainty of this estimate of the WIMP mass, we wish 
to demonstrate that keeping the terms oc r(Q mm ) in Eq.(45) indeed reduces the systematic 
error. This is demonstrated in Fig. 2, which shows the ratio of the reconstructed to input 
WIMP mass for experiments with "infinite" statistics. Here we have chosen X = 28 Si and 
Y = 76 Ge, and we have set Qmax = oo in order to isolate the systematic effect due to imperfect 
matching of the Q m in values in this figure.* Here the black (red) curves have been obtained with 
Qmin,Ge = Qmm, si = 10 (3) keV. The solid lines include the terms oc r(Q m [ n ), while the dashed 
ones do not. Clearly including these terms is beneficial: for m x > 20 GeV, the systematic shift 
with these terms included and Qmm = 10 keV is smaller than that without those terms with 
the much smaller Q mm = 3 keV. We also see that the systematic effect becomes large at small 
m x . This is not surprising, since a small m x implies a small Q ma x,kin, so that the lower cut on Q 
removes a correspondingly larger fraction of the total spectrum. The systematic error vanishes 
for m x = ^/mxTriY = 43 GeV, since for this value of m x we have ax = Q-y, i-e., i>i,x = i>i,y if 
Qmin,x = Qmin,Y- For larger m x the systematic effect changes sign, i.e., one now underestimates 
the true value of m x . However, if the terms oc r(Q min ) are included, the shift remains relatively 
small even for Qmm = 10 keV. If values Q min ^ 3 keV can be achieved, which should be possible 
[17], the systematic effect due to imperfect Q m in matching will be a concern only for very light 
WIMPs, m x < 20 GeV. 

Using Eq.(45) in Eq.(34) also leads to a lengthier expression for the statistical error on our 
estimate of m x . The first equation in (36) still holds, but the r(Q min ) terms give additional 
contributions to the final expression: 
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(46) 



Here we have introduced a short-hand notation for the six quantities on which the estimate of 
m x depends: 

Cl,X = In,X ; C 2 ,X = Io,X ; C 3j x = r X (Qmm,x) , (47) 

and similarly for the c^y. Estimators for cov(cj,Cj) have been given in Eqs.(26) and (31). 
Explicit expressions for the derivatives of TZ n with respect to these six quantities are collected 
in the Appendix. Note that a factor lZ n can be factored out of all these derivatives. With this 
factor moved out of the square brackets, the prefactor in Eq.(46) is identical to that in Eq.(36), 
i.e., its expectation value will again be given by k of Eq.(43). Of course, Eq.(46) reduces to 
Eq.(36) in the limit Q m m,x , Qmm,Y — ► 0. Note finally that Eq.(46) also holds for n = — 1, if 
the derivatives with respect to c±^x,y) are neglected, since in this case l n) x and l n y do not 
contribute to 1Z n given in Eq.(45). 

Finite lower energy cuts Q m m,(x,Y) can also easily be incorporated in the quantity TZ a ap- 
pearing in Eq.(40): 
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tMore exactly, for m x > 100 GeV we have set Qmax, Ge = 250 keV, in order to avoid complications due to the 
zero of the form factor of 76 Ge which occurs at Q ~ 270 keV. We have then matched Q max , si, so that the only 
systematic deviation of the reconstructed WIMP mass is indeed due to imperfect Q m i n matching. 



Correspondingly, Eq.(42) changes to 
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(49) 



where we have again used the short-hand notation of Eq.(47); note that c^x,y) does not appear 
here. Expressions for the derivatives of 1Z a are also given in the Appendix. 



3.3 Combined fit 



In the next step we wish to combine our estimators (34) for different n with each other, and with 
the estimator (40). This could be done via an overall covariance matrix describing the errors of 
these estimators and their correlations. The diagonal entries of this covariance matrix are given 
by Eqs.(46) and (49); the off-diagonal entries can be computed analogously. This would yield 
the overall best-fit value of m x as well as its Gaussian error. 

Here we pursue a slightly different procedure, based on a x 2 fit. This will yield the same 
best-fit value of m x , which we denote by m x>rec , but it has two advantages. First, X 2 ( m x,rcc) can 
be used as a measure of the quality of the fit, which in turn can be used to match the Q max values 
of the two experiments at least approximately. Secondly, it allows to determine asymmetric error 
intervals. Fig. 1 implies that the errors should indeed be asymmetric. For example, if the true 
m x is (much) larger than the masses of both target nuclei, the experimental upper bound one 
can derive will be quite large or even infinite, but one should still get a meaningful lower bound; 
the opposite is true if m x lies well below the mass of both target nuclei. 

We begin by defining fit functions 
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we analogously define n max + 2 functions f^y- Here n max determines the highest (generalized) 
moment of fi that is included in the fit. The fi are normalized such that they are dimensionless 
and very roughly of order unity; this alleviates numerical problems associated with the inversion 
of their covariance matrix. The first n max + 1 functions fi are basically our estimators (23) of the 
generalized moments defined in Eq.(12) with the term oc r(Q max ) omitted; as discussed earlier, 
the error on this quantity is likely to be so large that including this term will not be helpful. 
The last function is essentially the ratio appearing in Eq.(39). It is important to note that m x 
in Eqs.(50a) and (50b) is a fit parameter, not the true (input) value of the WIMP mass. Recall 
also that our estimator (25) for the integrals I n appearing in Eqs. (50a) and (50b) is independent 
of m x . Hence the first n max + 1 fit functions depend on m x only through the overall factor a 1 . 
The fi allow us to introduce a x 2 function: 
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Here C is the total covariance matrix. Since the X and Y quantities are statistically completely 
independent, C can be written as a sum of two terms: 

dj = cov (f iiX , fj,x) + cov (/ i>y , f jiY ) . (52) 

The entries of this matrix involving only the moments of the WIMP velocity distribution can 
be read off Eq.(82) of Ref. [14], with an obvious modification due to the normalization factor in 
Eq.(50a). Since the last fa can be computed from the same basic quantities, i.e., the counting 
rates at Q mm and the integrals Iq, the entries of the covariance matrix involving this last fit 
function can also be computed straightforwardly, using Eqs.(26)-(31). Of course, Eq.(51) can 
also be used to compute asymmetric error intervals from a single moment, by restricting the 
sum to a single term. 



4 Numerical results 

We are now ready to present some numerical results for the reconstructed WIMP mass. These 
results are based on Monte Carlo simulations of direct detection experiments. We assume that 
the scattering cross section is dominated by spin-independent interactions, and use the Woods- 
Saxon form for the elastic form factors F(Q) [16]. We describe the WIMP velocity distribution by 
a sum of a shifted Gaussian contribution [6] and a "late infall" component [12]; as in Ref. [14], for 
simplicity we model the latter as a 5— function, keeping the normalization N\± of this component 
as free parameter: 

h{v) = i^^l ( JL.) [e-e—M _ e -(^)>o 2 l + Nu s (v _ Vcsc) . (53 ) 



We take v = 220 km/s, v e = 1.05 1> § , and t> esc = 700 km/s. 

In Fig. 3 we show upper and lower bounds on the reconstructed WIMP mass, calculated from 
the requirement that x 2 exceeds its minimum by 1, for the idealized scenario where no cuts on Q 
have been applied. We took our default set of parameters, i.e., vanishing late infall component, 
and target nuclei X = 28 Si, Y = 76 Ge. This figure is based on simulating 2 x 5,000 experiments, 
where each experiment contains an expected 50 events; the actual number of events is Poisson- 
distributed around this expectation value. As mentioned earlier, taking equal numbers of events 
in both experiments minimizes the statistical error for fixed total number of events. As discussed 
in Ref. [14] , the error on the (high) moments is not quite Gaussian; the deviation becomes larger 
for smaller samples. The reason is that the high moments receive large contributions from the 
region of high Q, where on average very few events will lie; even the region where on average only 
a fraction of an event lies can contribute significantly. As a result, most simulated experiments 
will underestimate these moments, while a few (rare) experiments will overestimate them. In 
order to alleviate this problem, we only include moments up to n max = 2 in our fit. Moreover, 
we always show median, rather than mean, values for the (bounds on the) reconstructed WIMP 
mass. 

We see that the minus-first moment by itself leads to very poor bounds on m x . This is 
not surprising, since its error is dominated by the error on the counting rate at Qmin, which is 
determined only from the events in the first bin. The higher moments lead to increasingly tighter 
bounds. However, the higher moments are very strongly correlated. Also, the systematic effect 
due to the limited event samples discussed in the previous paragraph becomes larger for larger n; 



§ Strictly speaking, v e should oscillate around 1.05 vo with a period of one year [6]; we ignore this time 
dependence here. 




Figure 3: Median values of "1 <r" upper and lower bounds on the reconstructed WIMP mass in 
2 x 5, 000 simulated experiments with Si and Ge targets, as a function of the true value of m x . 
We generated on average 50 events per experiment. The short-dashed (black), dotted (green) 
and long-dashed (violet) curves show the upper and lower bounds on m x as determined from 
moments with n — —1,1 and 2, respectively; the dash-dotted (blue) curves labeled a show 
the corresponding limits derived from our assumption of equal cross sections for scattering on 
protons and neutrons. Finally, the solid (red) curves show the upper and lower bound as well as 
the median reconstructed m x for the global fit based on minimizing \ 2 of Eq.(51) with n max = 2. 
These results are for our standard halo, i.e., no late infall component, without any cuts on the 
recoil energy. 



for example, for n = 2 and a true m x = 1 TeV, the median upper end of the (nominal) 1 a range 
lies well below 1 TeV. The lower bound on m x derived from the assumption of equal scattering 
cross sections on protons and neutrons is similar to that derived from the first moment, but the 
corresponding upper bound is significantly worse. Nevertheless, this estimator of m x helps in 
narrowing down the error of the total fit, described by the upper and lower solid (red) curves. 
As expected from Fig. 1 the relative error on m x is minimal for m x = y/mxTny, although the 
increase towards smaller m x is less than expected from the behavior of the kinematical factor k 
alone. 

Unfortunately we also see that the median reconstructed m x starts to deviate from the 
input value if m x ^ 80 GeV. This is a direct consequence of the fact that the median value of 
the estimators of the higher moments is too small, as discussed above. For very large m x the 
median reconstructed WIMP mass even becomes independent of its true value; this is true also 
for the upper end of the error band. This systematic shift presents another argument in favor of 
imposing an upper cut Q max on Q, chosen sufficiently low that an average experiment will still 
have a few events not too far below Q max . 
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Figure 4: As in Fig. 3, except that we have imposed the cut Q max = 50 keV in both experiments. 
Note that the average of 50 events per experiment refers to the entire Q range, i.e., the number 
of events after cuts is smaller if Q max ,kin > Qmax- 



50 + 50 events, Si and Ge, standard halo, optimally matched O < 50 keV 




Figure 5: As in Fig. 4, except that we have imposed the cut Q max = 50 keV for the Ge target. 
The value of Q max for the Si target has been chosen such that it corresponds to the same WIMP 
velocity. 



For the case at hand, with rather small event samples, this would imply Q max ~ 50 keV. 
Unfortunately Fig. 4 shows that simply imposing the same cut on Q max in both targets makes 
the situation significantly worse. Now the median reconstructed m x starts to undershoot the 
input value already at m x ~ 60 GeV, and the asymptotic reconstructed m x for large input 
WIMP mass lies well below 100 GeV. Clearly we need to choose different Q max values for the 
two experiments if we want to get reliable results also for larger WIMP masses. 

Fig. 5 indicates that this should be possible, at least in principle. Here we have again applied 
a fixed upper cut Q max = 50 keV for the Ge experiment, but matched the cut on Q max for the 
Si experiment such that it corresponds to the same WIMP velocity: 

Qmax, Si ( I Qmax, Ge , (54) 

where a has been given in Eq.(5). We see that now the median reconstructed m x indeed tracks 
the input value even for very large WIMP masses. It should be noted that these results still only 
use 50 events on average per experiment before cuts. For example, for m x = 500 GeV Eq.(54) 
with Q maXj Ge = 50 keV implies Q maXj si = 21.7 keV, meaning that the expected number of events 
in the Si experiments is only about 21. This explains why the error bands at large values of m x 
are significantly wider here than in Fig. 3 where no cuts on the recoil energy have been imposed. 
Of course, the systematic deviation of the reconstructed m x from its true value makes the error 
bands in Fig. 3 largely meaningless for m x ^ 100 GeV. 

While the results in Fig. 5 look quite impressive, they suffer from the fact that optimal Q max 
matching, as in Eq.(54), is only possible if m x is already known. The left frame of Fig. 6 shows 
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Figure 6: Sensitivity of the \ 2 fit to the value of m x that is used as input in the Q max matching 
condition (54), assuming a true WIMP mass of 100 GeV, for our default choice of targets and 
fi{v). The solid (black) lines are for 2 x 50 events and Q max = 50 keV, while the dashed (red) 
curves assume 2 x 500 events and Q max = 100 keV; here Q max stands for the bigger of the 
values used for the two targets, the second one being fixed by Eq.(54). The curves terminate 
at m x ; 1VL = 25 GeV since for even smaller WIMP masses, Q m ax,kin < 50 keV, making matching 
superfluous. The left frame shows the median reconstructed WIMP mass from a \ 2 fit with 
n mffl = 2, and the right frame shows the corresponding mean value of X 2 (m XjI . ec ). 



that inputting the wrong value of m x into Eq.(54) will usually lead to a reconstructed WIMP 
mass in between this input value and the true value. Note that the median function of 

m Xi i n has a slope less than unity. This means that an iteration, where one starts with some input 
value of m x and uses the corresponding m x ^ ec as new input value into Eq.(54), will converge 
"on average". Unfortunately our Monte Carlo simulations show that for any one experiment, 
this procedure does not necessarily converge to a well-defined m Xirec ; rather, one often ends up 
in an endless loop over several values of m xrec . 

An alternative is to try an algorithm for Q max matching which is based on the minimum 
value of x 2 obtained in the fit; recall that this minimum defines the reconstructed WIMP mass. 
Unfortunately the right frame in Fig. 6 indicates that this may be difficult, at least with the 
initial small statistics expected. This figure shows that the mean value of x 2 ( m x,rcc) is almost 
independent of the value of m x input into Eq.(54) if one has only 2 x 50 events in the sample. 
The situation looks considerably more promising with 2 x 500 events, and larger allowed Q max : 
in this case X 2 ( m x,rec) has a clear, if rather broad, minimum where m Xtin equals the true WIMP 
mass, taken to be 100 GeV in this example. 

Note that this minimum has mean x 2 (m XiICC )/d.o.f. well below unity; the same is true for the 
entire solid curve. Here the number of degrees of freedom is three, since we fit four quantities 
with one free parameter. This indicates that our numerical procedure on average over-estimates 
the true errors somewhat. In fact, following Ref. [14], we add the "error on the error" to 
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Figure 7: As in Fig. 4, except that we have imposed the fixed cut Q max = 50 keV only on the 
Ge experiment and determined Q m ax,si by minimizing x 2 (m XjI . ec ). 



the diagonal entries of cov(J n ,J m ), in order to tame non-Gaussian tails in the distribution of 
measured moments of j\. 

One possibility is to choose the Q max values such that x 2 { m x,^c) is minimal. Note that 
this implies a double minimization: for fixed values of Q max ,x and Q max ,Y, ^x.rec is the WIMP 
mass that minimizes x 2 (m x ). In an outer loop, Qmax is varied. We found that varying both 
Qmax values leads to quite misleading results, especially for larger (true) WIMP masses and/or 
limited statistics. On the other hand, Fig. 4 shows that for small m x , Q mSuX matching is not 
really necessary; if m x > y/mxiTiY, the matching condition (54) implies that Q max for the lighter 
target should be smaller than that for the heavier target. 

In Fig. 7 we have therefore minimized X 2 { m x,rec) only for the lighter target, here Silicon. 
Note that we always require Q m ax,si < Qmax.Ge in this procedure; in Fig. 7 the latter is fixed to 
50 keV. We see that this leads to a systematic over-estimation of m x for small WIMP masses. 
For heavier WIMPs the results are somewhat better than for the case where both Qmax values 
are simply taken to be equal, see Fig. 4. However, this algorithm clearly still fails quite badly 
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Figure 8: Results for the median value of the reconstructed WIMP mass as well as the ends 
of its error interval. All results are based on the combined fit, using Eq.(51) with n max = 2. 
We assume Si and Ge targets with on average 50 events each before cuts, and took our default 
ansatz for f x . The dotted (green) lines show results for Q max ,si = Qmax,Ge = 100 keV, whereas 
the solid (black) lines have been obtained using Eq.(54) with the bigger of the two Q m ax values 
fixed to 100 keV. Finally, the dashed (red) lines are for the case that Q ma x,Ge = 100 keV, whereas 
Qmax,Si has been chosen such that x 2 ( m x,rec) is minimal. 



for m x > 100 GeV. This is partly due to the low maximal value of Qmax assumed here. The fact 
that the kinematic upper bound Q max ,kin of Eq. (8) increases with m x indicates that the region of 
large Q is more sensitive to the value of m x if the WIMP mass exceeds the mass of the heaviest 
target nucleus. 

In Fig. 8 we have therefore increased the cut Qmax to 100 keV; of course, this cut only becomes 
effective once Q max ,kin > 100 keV. Unlike in the previous figures, we here show results only for 
the final fit; the values of m x derived from single observables are no longer shown. This allows 
us to show results for three different choices of <5max,Si an d Q m ax,Ge- The dotted (green) curves 
show the median reconstructed WIMP mass and its "lcr" upper and lower bounds for the case 
where both Q max values have been fixed to 100 keV. Due to the higher Q max chosen here, this 
works for considerably larger WIMP masses than in the corresponding Fig. 4, where both Qmax 
values had been fixed to 50 keV. However, for large m x the median reconstructed WIMP mass 
is still only slightly above 100 GeV. The solid (black) lines show results for the case that perfect 
Qmax matching has been applied, using Eq.(54). Comparison with Fig. 5 shows that increasing 
Qmax actually slightly increases the width of the errors on the reconstructed WIMP mass. The 
reason is that Q max is now so large that the median values of the estimators for the (generalized) 
moments of f\ fall somewhat below the true values. Nevertheless the results for this optimal 
Qmax matching remain very encouraging. 

Most importantly, our simple algorithm of fixing Qmax,Ge, in this case to 100 keV, and de- 
termining <5max,Si by minimizing x 2 (m XiI . cc ) now also seems to work reasonably well for WIMP 
masses up to ~ 500 GeV. For m x ^ 100 GeV the median WIMP mass determined in this way 
again over-estimates its true value by 15 to 20%; however, the median value of the "la" lower 
bound lies below the true WIMP mass for all values of m x . Similarly, the median "lcr" upper 
bound now lies well above the true WIMP mass even for m x = 1 TeV, in sharp contrast to the 
results shown in Fig 7. 

We had argued earlier, based on the results of Fig. 1, that Si and Ge is likely to be close to 
the optimal choice among the target nuclei currently being employed. In Fig. 9 we back this up 
by showing results analogous to those of Fig. 8, but with Argon and Xenon targets. Since the 
elastic form factor of 136 Xe has a zero near 95 keV, we have lowered the upper bound on Q ma _ x 
to 75 keV; larger values would start to probe the very sparsely populated region near the zero, 
leading to too small median values of the estimated moments, whereas smaller values would lead 
to even worse sensitivity to large WIMP masses. We see that, except for small WIMP masses, 
the results are clearly worse than those shown in Fig. 8. 

The starting point of our discussion was that we wanted to devise ways to estimate the WIMP 
mass which are independent of any assumptions on the WIMP velocity distribution f x . In order 
to test this, we have simulated Si and Ge experiments, allowing 25% of the local WIMP flux to 
come from a "late infall" component, i.e., we fixed N\± = 0.25 in Eq.(53). The results are shown 
in Fig. 10. We see that the results are very similar to those with N\± = shown in Fig. 8 if optimal 
Qmax matching (54) is used. On the other hand, the results for non-optimal choices of Q max get 
somewhat worse. This is true both for the simple choice Q maXi si = Qmax,Go where significant 
deviations set in for lower values of the WIMP mass, and for our algorithm of determining 
Qmax,si by minimizing x 2 { m x,rec)- I n the latter case, the systematic difference between median 
reconstructed and true WIMP mass is larger than for N\± = 0, both for small and for large m x . 
The reason for this degradation is that introducing a large late-infall component, corresponding 
to WIMPs with velocity about three times larger than the mean velocity of the shifted Gaussian 
component, increases the number of events at large recoil energy. Hence correct Q max matching 
becomes more important. Note, however, that the true m x always lies within the median limits 
of the "la" error interval estimated from our algorithmic Q max matching. 
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Figure 9: As in Fig. 8, except for 40 Ar and 136 Xe targets, and with both Q max values restricted 
to not exceed 75 keV. 



50 + 50 events, Si and Ge, halo with 25% late infall, Q < 100 keV 
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Figure 10: As in Fig. 8, except that we allowed a 25% late infall component in the local WIMP 
flux, by setting N u . = 0.25 in Eq.(53). 




Figure 11: As in Fig. 8, except for 2 x 500 events before cuts. 



So far we have assumed that each experiment "only" has an exposure corresponding to 50 
events before cuts. In Fig. 11 we raise this number by a factor of 10. Not surprisingly, all error 
bars shrink by a factor ^ 3 compared to the situation of Fig. 8. The error interval also remains 
approximately symmetric out to much larger values of m x . However, the larger number of events 
does not significantly change the median reconstructed m x if one simply fixes both Q max values 
to 100 keV. This is not surprising: for large m x this implies that V2,si and f 2,Ge in the definition 
(12) of the generalized moments, and hence the moments themselves, are quite different. Hence 
the estimators of these moments will not agree if the true WIMP mass is used. On the other 
hand, our algorithm for fixing Q max ,Ge now seems to perform very well over the entire range of 
WIMP masses shown. In particular, the median reconstructed WIMP mass now overshoots its 
true value by only a few percent for m x ^ 100 GeV, and remains close to the true value even at 
m x = 1 TeV. Unfortunately we will see shortly that for large WIMP masses our algorithm still 
has some problems even in this case. 

The problem lies in the distribution of the reconstructed WIMP masses in the simulated 
experiments. This distribution is supposed to be characterized by the error intervals shown in 



Figs. 3-5 and 7—11. In order to see how well this works, we introduce the quantity 
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Here m x is the true WIMP mass, m 



its reconstructed value, m Xi i i( 2 ) is the "1 (2) a" lower 
bound satisfying X 2 ( m xM 1 2)) = X 2 ( m x,rec) + 1 (4), and m Xim i(2) are the corresponding "1 (2) cr" 
upper bounds. In the limit of purely Gaussian errors, where x 2 °f Eq.(51) is simply a parabola, 
{8m) 2 would itself be a x 2 variable, measuring the difference between the true and the recon- 
structed WIMP mass in units of the error of the reconstruction. However, we saw earlier that 
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Figure 12: Normalized distribution of the variable 8m defined in Eq.(55) for 5,000 simulated 
experiments, for true WIMP mass m x = 50 GeV. The other parameters, are as in Fig. 8. The 
solid (black) histogram shows results for perfect Qmax matching with Q max < 100 keV, based on 
Eq.(54), whereas the dotted (green) histogram is for fixed equal Q max values of 100 keV. The 
dashed (red) histogram shows results when Q max ,si is determined by minimizing x 2 ( m x,rec)- 



the error intervals are often quite asymmetric. Similarly, the distance between the "2cr" and 
"la" limits can be quite different from the distance between the "la" limit and the central 
value. The definition (55) takes these differences into account, and also keeps track of the sign 
of the deviation: if the reconstructed WIMP mass is larger (smaller) than the true one, 5m is 
positive (negative). Moreover, \8m\ < 1 (2) if and only if the true WIMP mass lies between the 
"experimental" 1 (2) a limits. 

In Fig. 12 we show the distribution of 5m calculated from 5,000 simulated experiments, 
assuming a rather small WIMP mass, m x = 50 GeV. The other parameters have been fixed 
as in Fig. 8. In this case simply fixing both Q max values to 100 keV still works fine, since the 
kinematic maximum values of Q lie only slightly above 100 keV (at 122 keV for Si and 132 
keV for Ge). The distributions for fixed Q max , or for optimal Q max matching, look somewhat 
lopsided, since the error interval is already asymmetric, with m xh n — ^ x ,rec > m x,rec — m x ,ioi- 
As a result, negative values of 5m have a larger denominator than positive values, hence the 
distribution is narrower for 5m < 0. These distributions also indicate that our errors are indeed 
over-estimated, since nearly 90% of the simulated experiments have \5m\ < 1; we remind the 
reader that a usual la interval should only contain some 68% of the experiments. 

We saw in Fig. 8 that our algorithm for determining Q max ,si tends to overestimate the WIMP 
mass if the latter is small. This is reflected by the dashed (red) histogram in Fig. 12, which has 
significantly more entries at positive values than at negative values. Moreover, this histogram 
is rather flat between 5m = 1 and 2. Since our algorithm is based on a double minimization 
of x 2 defined in Eq.(51), it is not very surprising that the resulting final x 2 (m Xirec ) values are 
distributed quite differently from what one finds after a single minimization step. Nevertheless 
it is reassuring that some 70% of the simulated experiments lead to \5m\ < 1. 
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Figure 13: Distribution of 5m defined in Eq.(55) calculated from 5,000 simulated experiments. 
Parameters and notations are as in Fig. 12, except that we have increased the WIMP mass to 
200 GeV. In the right frame we have in addition increased the average number of events (before 
cuts) in each experiment from 50 to 500. Note that the bins at 5m = ±5 are overflow bins, i.e., 
they also contain all experiments with \5m\ > 5. 



Unfortunately Fig. 13 shows that the situation becomes much less favorable at the larger 
WIMP mass of 200 GeV. Here we show results with 50 (left) and 500 (right) events per ex- 
periment (before cuts), with the other parameters as in Fig. 12. We see that for optimal Qmax 
matching a large majority of the simulated experiments still satisfy \5m\ < 1. The fact that this 
fraction decreases from nearly 90% for the smaller events samples to about 85% for the larger 
samples indicates that our error estimates become a bit more reliable for larger event numbers.^ 
Note also that the secondary peak at 5m ~ — 1 is due to the change of denominator in the 
definition (55). 

We already saw in Figs. 8 and 11 that setting Qmax.si = Qmax.Ge = 100 keV significantly 
under-estimates the true WIMP mass if the latter exceeds 100 GeV. This is borne out by 
Figs. 13. Since the statistical error decreases with increasing number of events, 5m is much 
smaller in the right frame; in fact, most of the simulated experiments now fall in the "underflow 
bin" 5m = —5, which also contains all experiments giving even smaller values. In other words, 
most of the simulated experiments would give a reconstructed WIMP mass more than five 
estimated standard deviations below the true value, if no Q max matching is used. 

Unfortunately the error estimates resulting from our algorithm of determining Q max ,si by 
minimizing x 2 { m x,rcc) are also not very reliable in this case. For the smaller event sample, we 
find that about 58% of the simulated experiments yield \5m\ < 1, while nearly all experiments 
give 1 5m \ < 2. While these numbers are not so different from the corresponding Gaussian 
predictions, the distribution of 5m is clearly highly non-Gaussian in this case. Just as in the 
scenario without Q max matching the error estimates actually become less reliable with increasing 
event samples. If both experiments contain on average 500 events each, less than 40% of the 
experiments have \5m\ < 1, while more than 30% of the simulations yield |<5m| > 2. In fact, 
the most likely value of 5m is now close to 3. In view of these observations, the fact that the 
median 5m is close to zero, so that the median reconstructed WIMP mass is close to the true 
value as already shown in Fig. 11, seems almost accidental. Recall also that the nominal "lcr" 
uncertainty of the reconstructed WIMP mass still amounts to about 40 GeV in this case. This 
means that the most likely value of m XiVec predicted by our algorithm exceeds the true value by 
more than 100 GeV. This clearly leaves some room for improvements. The fact that optimal 
Qmax matching continues to give good results for both the reconstructed WIMP mass and its 
error indicates that better data-based algorithms might very well exist. 



5 Summary and Conclusions 

In this paper we described methods to determine the mass of a Weakly Interacting Massive 
Particle detected "directly" , i.e., through the recoil energy deposited in a detector by the recoiling 
nucleus after a WIMP scattered elastically off this nucleus. Our methods are model-independent 
in the sense that they do not need any assumption about the WIMP velocity distribution. The 
price one has to pay for this is that one will need positive signals in at least two different 
detectors, employing different target nuclei. 

Our methods are based on our earlier work [14] on reconstructing the WIMP velocity dis- 
tribution, which we briefly reviewed in Sec. 2. In this earlier work, which was based on results 
from a single (simulated) experiment, the WIMP mass m x was an input. In Sec. 3 we showed 
how one can determine m x by equating results obtained by different experiments. Here the 
moments of the velocity distribution function are particularly useful, since all events in the sam- 

" Assuming, unrealistically, that there are 50,000 events in each experiment, we find that only about 75% of 
all experiments have \Sm\ < 1, indicating a very slow approach to the Gaussian value of about 68%. 



pie contribute to any given moment, leading to relatively low statistical uncertainties. We also 
described a method for determining m x that can be used if the ratio of WIMP scattering cross 
sections on protons and neutrons is known; this is true, for example, for the spin-independent 
scattering of supersymmetric neutralinos, where these two cross sections are nearly equal. We 
also showed how to combine these methods using a x 2 fitting procedure. 

Sec. 4 was devoted to a detailed numerical analysis of our methods. We saw that, assuming 
the sizes of the event samples are fixed, the statistical errors will be smaller for larger mass 
difference between the two target nuclei. In practice experiments with heavier targets will 
accumulate more events, assuming equal exposure, at least if the spin-independent contribution 
to the scattering cross section dominates. However, the number of useful events (after cuts, 
preferably in an almost background-free energy range) also depends on other factors, besides 
the masses of the target nuclei. 

In our discussion we saw that, for WIMP masses exceeding 100 GeV or so, the maximal recoil 
energy Qmax of accepted signal events plays a crucial role. Existing experiments have Q max < 100 
keV. If both targets used fixed Q max values of this order or even smaller, a significant systematic 
error on the extracted WIMP mass results. In principle this problem can be solved by matching 
the Qmax values of the two experiments. The problem is that perfect matching requires prior 
knowledge of the WIMP mass. We tried two algorithms to overcome this problem. Determining 
Qmax of one experiment iteratively should converge "on average", but in a given experiment 
often leads to an endless loop, rather than a specific value of Q max ; this problem is particularly 
severe for small event samples. On the other hand, determining Q max of the experiment with 
the lighter target nucleus by minimizing x 2 also with respect to this quantity over-estimates the 
WIMP mass if it is small, and leads to unreliable error estimates if the WIMP mass is larger, the 
problem becoming worse with increasing event samples. However, the fact that optimal Q max 
matching works well in all cases, for both the median reconstructed WIMP mass and its error 
(which tends to be over-estimated by our expressions), gives us hope that a better algorithm for 
Qmax matching can be found which only relies on the data. One possibility that might be worth 
exploring is to employ a combination of an iterative procedure and a second x 2 minimization, 
where the latter is used only if the former does not converge to a well-defined value of the 
reconstructed WIMP mass. 

We also found that imposing a cut Q ma , x may actually be beneficial for small event samples. 
This is related to our earlier observation [14] that a typical experiment will under-estimate the 
higher moments of f\, which receive significant contributions from recoil energies where only a 
fraction of an event is expected to occur in a given experiment. This problem becomes more 
acute for heavier target nuclei, since they have softer form factors. In particular, using Xenon 
rather than Germanium does not improve the determination of m x , since the spin-independent 
elastic form factor of Xenon as predicted by the Woods-Saxon ansatz vanishes for Q ~ 95 keV. 
In contrast, the lower (threshold) energy of the experiment does not seem to be very important, 
if it can be pushed down to values near 3 keV or less. 

Our analysis is idealized in that we ignore backgrounds, systematic uncertainties as well 
as the finite energy resolution. The relative error on the recoil energy in existing experiments 
is small compared to our most optimistic relative WIMP mass error estimates even with 500 
events per experiment, so ignoring it should be a good approximation. Modern methods of 
discriminating between nuclear recoils and other events, combined with muon veto and good 
shielding, hold out the possibility of keeping (some) experiments nearly background-free also in 
future. 

In our numerical analysis we have ignored the expected annual modulation of the WIMP 
flux. In practice this can be done if one simply sums all events over (at least) one full calendar 



year. In principle one can also use our methods for subsets of data collected during specific 
times of the year. However, at least if the standard "shifted Gaussian" velocity distribution is 
approximately correct, we do not expect the small annual modulation to play a significant role, 
even if one compares experiments taken during different parts of the year. 

We saw that our methods work best if the WIMP mass lies in between the masses of the two 
target nuclei. Even in that case the error will likely be significantly larger than the error on m x 
from collider experiments, if the WIMP is part of a well-motivated extension of the Standard 
Model of particle physics, e.g., if it is the lightest neutralino [6] or the lightest T— odd particle 
in "Little Higgs" models [18]. It will nevertheless be crucial to determine the WIMP mass from 
direct and/or indirect Dark Matter experiments as precisely as possible, in order to make sure 
that the particle produced at colliders is indeed the WIMP detected by these experiments. Once 
one is confident of this identification, one can use further collider measurements to constrain 
the WIMP couplings. This in turn will allow to calculate the WIMP-nucleus scattering cross 
section. Together with the determination of the WIMP velocity distribution [14], this will 
then yield a determination of the local WIMP number density via the total counting rate in 
direct detection experiments. Knowledge of the WIMP couplings will also permit prediction 
of the WIMP annihilation cross section. Together with the Dark Matter density inferred from 
cosmological observations, this will allow to test our understanding of the early universe [19]. 
A determination of the WIMP mass from Dark Matter detection experiments is thus a crucial 
ingredient in many analyses that shed light on the dark sector of the universe. 
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A Derivatives needed in the error analysis 



At the end of Sec. 2 we gave the covariance matrix of the quantities appearing in the definition 
of the lZ n as well as 1Z a . In Eqs.(46) and (49) we also gave expressions relating the errors on the 
reconstructed WIMP masses to the errors on 1Z n and 1Z a . The only missing ingredients in the 
calculation of the errors on our various estimators of m x are the first derivatives of lZ n and lZ a . 
We begin with the former. From Eq.(45), it can be found directly that 
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By first exchanging Q^^ x 2 an d (n+l)I n ,x with Q l J^ nX an d Io,x, respectively, and then replacing 
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Note that a factor 7£ n appears in all these expressions; this allows to cast the final result for 
the error on m x estimated using moments of f\ into a form analogous to that in Eq.(36) even 
in the presence of non-trivial cuts on the recoil energy, with the same prefactor. Moreover, 
all the 7 ,x, 7 y, I n ,x, I n ,Y should be understood to be computed according to Eq.(14) or its 
discretization (25) with integration limits Qmin and Qmax specific for that target. 
Similarly, the derivatives of lZ a can be computed from Eq.(48): 
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The derivatives with respect to the Y variables can be obtained from Eqs.(A3a) and (A3b) by 
simply changing X — > Y everywhere and changing the overall plus signs to minus signs. 
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